home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / quadpack / dqelg.f < prev    next >
Text File  |  1996-07-19  |  6KB  |  185 lines
  1.       SUBROUTINE DQELG(N,EPSTAB,RESULT,ABSERR,RES3LA,NRES)
  2. C***BEGIN PROLOGUE  DQELG
  3. C***REFER TO  DQAGIE,DQAGOE,DQAGPE,DQAGSE
  4. C***ROUTINES CALLED  D1MACH
  5. C***REVISION DATE  830518   (YYMMDD)
  6. C***KEYWORDS  EPSILON ALGORITHM, CONVERGENCE ACCELERATION,
  7. C             EXTRAPOLATION
  8. C***AUTHOR  PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
  9. C           DE DONCKER,ELISE,APPL. MATH & PROGR. DIV. - K.U.LEUVEN
  10. C***PURPOSE  THE ROUTINE DETERMINES THE LIMIT OF A GIVEN SEQUENCE OF
  11. C            APPROXIMATIONS, BY MEANS OF THE EPSILON ALGORITHM OF
  12. C            P.WYNN. AN ESTIMATE OF THE ABSOLUTE ERROR IS ALSO GIVEN.
  13. C            THE CONDENSED EPSILON TABLE IS COMPUTED. ONLY THOSE
  14. C            ELEMENTS NEEDED FOR THE COMPUTATION OF THE NEXT DIAGONAL
  15. C            ARE PRESERVED.
  16. C***DESCRIPTION
  17. C
  18. C           EPSILON ALGORITHM
  19. C           STANDARD FORTRAN SUBROUTINE
  20. C           DOUBLE PRECISION VERSION
  21. C
  22. C           PARAMETERS
  23. C              N      - INTEGER
  24. C                       EPSTAB(N) CONTAINS THE NEW ELEMENT IN THE
  25. C                       FIRST COLUMN OF THE EPSILON TABLE.
  26. C
  27. C              EPSTAB - DOUBLE PRECISION
  28. C                       VECTOR OF DIMENSION 52 CONTAINING THE ELEMENTS
  29. C                       OF THE TWO LOWER DIAGONALS OF THE TRIANGULAR
  30. C                       EPSILON TABLE. THE ELEMENTS ARE NUMBERED
  31. C                       STARTING AT THE RIGHT-HAND CORNER OF THE
  32. C                       TRIANGLE.
  33. C
  34. C              RESULT - DOUBLE PRECISION
  35. C                       RESULTING APPROXIMATION TO THE INTEGRAL
  36. C
  37. C              ABSERR - DOUBLE PRECISION
  38. C                       ESTIMATE OF THE ABSOLUTE ERROR COMPUTED FROM
  39. C                       RESULT AND THE 3 PREVIOUS RESULTS
  40. C
  41. C              RES3LA - DOUBLE PRECISION
  42. C                       VECTOR OF DIMENSION 3 CONTAINING THE LAST 3
  43. C                       RESULTS
  44. C
  45. C              NRES   - INTEGER
  46. C                       NUMBER OF CALLS TO THE ROUTINE
  47. C                       (SHOULD BE ZERO AT FIRST CALL)
  48. C
  49. C***END PROLOGUE  DQELG
  50. C
  51.       DOUBLE PRECISION ABSERR,DABS,DELTA1,DELTA2,DELTA3,DMAX1,D1MACH,
  52.      *  EPMACH,EPSINF,EPSTAB,ERROR,ERR1,ERR2,ERR3,E0,E1,E1ABS,E2,E3,
  53.      *  OFLOW,RES,RESULT,RES3LA,SS,TOL1,TOL2,TOL3
  54.       INTEGER I,IB,IB2,IE,INDX,K1,K2,K3,LIMEXP,N,NEWELM,NRES,NUM
  55.       DIMENSION EPSTAB(52),RES3LA(3)
  56. C
  57. C           LIST OF MAJOR VARIABLES
  58. C           -----------------------
  59. C
  60. C           E0     - THE 4 ELEMENTS ON WHICH THE COMPUTATION OF A NEW
  61. C           E1       ELEMENT IN THE EPSILON TABLE IS BASED
  62. C           E2
  63. C           E3                 E0
  64. C                        E3    E1    NEW
  65. C                              E2
  66. C           NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW
  67. C                    DIAGONAL
  68. C           ERROR  - ERROR = ABS(E1-E0)+ABS(E2-E1)+ABS(NEW-E2)
  69. C           RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST VALUE
  70. C                    OF ERROR
  71. C
  72. C           MACHINE DEPENDENT CONSTANTS
  73. C           ---------------------------
  74. C
  75. C           EPMACH IS THE LARGEST RELATIVE SPACING.
  76. C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
  77. C           LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON
  78. C           TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER
  79. C           DIAGONAL OF THE EPSILON TABLE IS DELETED.
  80. C
  81. C***FIRST EXECUTABLE STATEMENT  DQELG
  82.       EPMACH = D1MACH(4)
  83.       OFLOW = D1MACH(2)
  84.       NRES = NRES+1
  85.       ABSERR = OFLOW
  86.       RESULT = EPSTAB(N)
  87.       IF(N.LT.3) GO TO 100
  88.       LIMEXP = 50
  89.       EPSTAB(N+2) = EPSTAB(N)
  90.       NEWELM = (N-1)/2
  91.       EPSTAB(N) = OFLOW
  92.       NUM = N
  93.       K1 = N
  94.       DO 40 I = 1,NEWELM
  95.         K2 = K1-1
  96.         K3 = K1-2
  97.         RES = EPSTAB(K1+2)
  98.         E0 = EPSTAB(K3)
  99.         E1 = EPSTAB(K2)
  100.         E2 = RES
  101.         E1ABS = DABS(E1)
  102.         DELTA2 = E2-E1
  103.         ERR2 = DABS(DELTA2)
  104.         TOL2 = DMAX1(DABS(E2),E1ABS)*EPMACH
  105.         DELTA3 = E1-E0
  106.         ERR3 = DABS(DELTA3)
  107.         TOL3 = DMAX1(E1ABS,DABS(E0))*EPMACH
  108.         IF(ERR2.GT.TOL2.OR.ERR3.GT.TOL3) GO TO 10
  109. C
  110. C           IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
  111. C           ACCURACY, CONVERGENCE IS ASSUMED.
  112. C           RESULT = E2
  113. C           ABSERR = ABS(E1-E0)+ABS(E2-E1)
  114. C
  115.         RESULT = RES
  116.         ABSERR = ERR2+ERR3
  117. C ***JUMP OUT OF DO-LOOP
  118.         GO TO 100
  119.    10   E3 = EPSTAB(K1)
  120.         EPSTAB(K1) = E1
  121.         DELTA1 = E1-E3
  122.         ERR1 = DABS(DELTA1)
  123.         TOL1 = DMAX1(E1ABS,DABS(E3))*EPMACH
  124. C
  125. C           IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT
  126. C           A PART OF THE TABLE BY ADJUSTING THE VALUE OF N
  127. C
  128.         IF(ERR1.LE.TOL1.OR.ERR2.LE.TOL2.OR.ERR3.LE.TOL3) GO TO 20
  129.         SS = 0.1D+01/DELTA1+0.1D+01/DELTA2-0.1D+01/DELTA3
  130.         EPSINF = DABS(SS*E1)
  131. C
  132. C           TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND
  133. C           EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE
  134. C           OF N.
  135. C
  136.         IF(EPSINF.GT.0.1D-03) GO TO 30
  137.    20   N = I+I-1
  138. C ***JUMP OUT OF DO-LOOP
  139.         GO TO 50
  140. C
  141. C           COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST
  142. C           THE VALUE OF RESULT.
  143. C
  144.    30   RES = E1+0.1D+01/SS
  145.         EPSTAB(K1) = RES
  146.         K1 = K1-2
  147.         ERROR = ERR2+DABS(RES-E2)+ERR3
  148.         IF(ERROR.GT.ABSERR) GO TO 40
  149.         ABSERR = ERROR
  150.         RESULT = RES
  151.    40 CONTINUE
  152. C
  153. C           SHIFT THE TABLE.
  154. C
  155.    50 IF(N.EQ.LIMEXP) N = 2*(LIMEXP/2)-1
  156.       IB = 1
  157.       IF((NUM/2)*2.EQ.NUM) IB = 2
  158.       IE = NEWELM+1
  159.       DO 60 I=1,IE
  160.         IB2 = IB+2
  161.         EPSTAB(IB) = EPSTAB(IB2)
  162.         IB = IB2
  163.    60 CONTINUE
  164.       IF(NUM.EQ.N) GO TO 80
  165.       INDX = NUM-N+1
  166.       DO 70 I = 1,N
  167.         EPSTAB(I)= EPSTAB(INDX)
  168.         INDX = INDX+1
  169.    70 CONTINUE
  170.    80 IF(NRES.GE.4) GO TO 90
  171.       RES3LA(NRES) = RESULT
  172.       ABSERR = OFLOW
  173.       GO TO 100
  174. C
  175. C           COMPUTE ERROR ESTIMATE
  176. C
  177.    90 ABSERR = DABS(RESULT-RES3LA(3))+DABS(RESULT-RES3LA(2))
  178.      *  +DABS(RESULT-RES3LA(1))
  179.       RES3LA(1) = RES3LA(2)
  180.       RES3LA(2) = RES3LA(3)
  181.       RES3LA(3) = RESULT
  182.   100 ABSERR = DMAX1(ABSERR,0.5D+01*EPMACH*DABS(RESULT))
  183.       RETURN
  184.       END
  185.